#include"../header/monteCarlo.h"

#define PPI 10000000
#define MAX_ITERACOES 200
#define N_THREADS 250

void monteCarloSeq(FILE *arquivo)
{	
	srand48(time(NULL));

	/*declaracao*/
	int j = 1;
	int bits = 33219281; /* 10^7 * log(10) */ 
	double x, y;
	mpz_t count, total;
	unsigned long int count_aux, i;	

	struct timeval inicial, final;

	mpf_t pi0, pi;
	mpf_t aux, delta, precision;

	/*inicializacao*/
	printf("inicializacao...\n");	
	mpf_set_default_prec(bits);
	
	mpf_init(pi);
	mpf_init(pi0);
	mpf_init(delta);
	mpf_init(precision);
	mpf_init(aux);

	mpz_init(count);
	mpz_init(total);

	/*setando a precisao*/
	printf("setando variaveis...\n");

	inicial = tempo();

	mpf_set_d(precision,10);
	mpf_pow_ui(precision,precision,10000000);
	mpf_ui_div(precision,1,precision);
	

	printf("iniciando procedimento...\n");

	do
	{	
		count_aux = 0;
		for(i = 0; i < PPI; ++i)
		{
			x = drand48();
			y = drand48();
			
			if( (x*x + y*y) <= 1.0f )
				count_aux++;			
		}

		mpz_add_ui(total,total,(unsigned long int) PPI);
		mpz_add_ui(count,count, count_aux);
		
		/*calculo do pi*/
		mpf_set_z(pi,count);
		mpf_set_z(aux,total);
		mpf_div(pi,pi,aux);
		mpf_mul_ui(pi,pi,(unsigned long int) 4);

		/*comparacao*/
		mpf_sub(delta,pi,pi0);
		mpf_abs(delta,delta);	

		/*atualizando pi0*/
		mpf_set(pi0,pi);
		printf("fim da iteracao %d\n",j);
		j++;

	} while((mpf_cmp(delta,precision) > 0) && j < MAX_ITERACOES);

	final = tempo();

    printf("Tempo de execucao: %.3lf s\n", difftempo(final,inicial));

	mpf_out_str (arquivo, 10, 10000001, pi);

	mpf_clear(pi);
	mpf_clear(pi0);
	mpf_clear(aux);
	mpf_clear(delta);
	mpf_clear(precision);

	mpz_clear(count);
	mpz_clear(total);

}

unsigned long int count_aux_par[N_THREADS];

int j_par;

void * calc_rand(void * n)
{
	int t = (int) n;
	struct drand48_data buffer;
	unsigned long int i;
	double x, y;

	srand48_r(time(NULL)+10000000*j_par+1000*t,&buffer);

	count_aux_par[(int)n] = 0;
	
	for(i = 0; i < PPI/(int)N_THREADS; ++i)
	{
		drand48_r(&buffer,&x);
		drand48_r(&buffer,&y);

		if( (x*x + y*y) <= 1.0f )
			count_aux_par[(int)n]++;			
	}

	return NULL;
}


void monteCarloPar(FILE *arquivo)
{

	/*declaracao*/
	int i;
	int bits = 33219281; /* 10^7 * log(10) */ 
	mpz_t count, total;

	pthread_t gerador[N_THREADS];

	struct timeval inicial, final;

	mpf_t pi0, pi;
	mpf_t aux, delta, precision;

	/*inicializacao*/
	printf("inicializacao...\n");	
	mpf_set_default_prec(bits);

	j_par = 1;
	
	mpf_init(pi);
	mpf_init(pi0);
	mpf_init(delta);
	mpf_init(precision);
	mpf_init(aux);

	mpz_init(count);
	mpz_init(total);

	/*setando a precisao*/
	printf("setando variaveis...\n");

	inicial = tempo();

	mpf_set_d(precision,10);
	mpf_pow_ui(precision,precision,10000000);
	mpf_ui_div(precision,1,precision);
	

	printf("iniciando procedimento...\n");

	do
	{
		for(i = 0; i < N_THREADS; i++)
			pthread_create(&(gerador[i]), NULL, calc_rand, (void *) i);


		for(i = 0; i < N_THREADS; i++)
			pthread_join(gerador[i], NULL);

		for(i = 0; i < N_THREADS; i++)
			mpz_add_ui(count,count, count_aux_par[i]);

		mpz_add_ui(total,total,(unsigned long int) PPI);
		mpf_set_z(aux,total);

		/*calculo do pi*/
		mpf_set_z(pi,count);
		mpf_div(pi,pi,aux);
		mpf_mul_ui(pi,pi,(unsigned long int) 4);

		/*comparacao*/
		mpf_sub(delta,pi,pi0);
		mpf_abs(delta,delta);	

		/*atualizando pi0*/
		mpf_set(pi0,pi);
		printf("fim da iteracao %d\n",j_par);
		j_par++;

	} while((mpf_cmp(delta,precision) > 0) && j_par < MAX_ITERACOES);

	final = tempo();

    printf("Tempo de execucao: %.3lf s\n", difftempo(final,inicial));

	mpf_out_str (arquivo, 10, 10000001, pi);

	mpf_clear(pi);
	mpf_clear(pi0);
	mpf_clear(aux);
	mpf_clear(delta);
	mpf_clear(precision);

	mpz_clear(count);
	mpz_clear(total);


}

